Genome-wide association analysis identify candidate genes for feed efficiency and growth traits in Wenchang chickens

Background Wenchang chickens are one of the most popular local chicken breeds in the Chinese chicken industry. However, the low feed efficiency is the main shortcoming of this breed. Therefore, there is a need to find a more precise breeding method to improve the feed efficiency of Wenchang chickens. In this study, we explored important candidate genes and variants for feed efficiency and growth traits through genome-wide association study (GWAS) analysis. Results Estimates of genomic heritability for growth and feed efficiency traits, including residual feed intake (RFI) of 0.05, average daily food intake (ADFI) of 0.21, average daily weight gain (ADG) of 0.24, body weight (BW) at 87, 95, 104, 113 days of age (BW87, BW95, BW104 and BW113) ranged from 0.30 to 0.44. Important candidate genes related to feed efficiency and growth traits were identified, such as PLCE1, LAP3, MED28, QDPR, LDB2 and SEL1L3 genes. Conclusion The results identified important candidate genes for feed efficiency and growth traits in Wenchang chickens and provide a theoretical basis for the development of new molecular breeding technology. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10559-w.


Introduction
Poultry production is an important entreprise worldwide.Chicken is considered a healthy white meat source that is lower in fat, calories and cholesterol than other red meat sources [1].In recent years, the promotion of white meat consumption has gradually become a trend [2].As an important part of the white meat market, chickens represent an efficient and inexpensive source of animal protein [3].Since the 1980s, with the continuous development of modern broiler production, China has become the second largest country in the world in terms of chicken production and consumption [4].The main breeds of chickens are native breeds and broilers [5].Wenchang chicken (Fig. 1) is one of the most popular chicken breeds in the local chicken industry, originating from Hainan Island in the South China Sea and has been raised for 400 years.It is famous for its excellent meat quality and is one of the four famous dishes from Hainan [6].It's an economic mainstay of animal husbandry in Hainan Province, with an annual output of nearly 100 million chickens and a total output value of 1.78 billion dollars in 2020 [7].
However, the primary shortcoming of Wenchang chickens is their low feed efficiency [8].Therefore, there is an urgent need to find a more reliable breeding method to improve the feed efficiency of Wenchang chickens.Feed represents more than 70% of the total cost of poultry production, and improving feed efficiency has consistently been the goal of any chicken breeding strategy [9].Feed efficiency is contingent upon the relation between the feed intake (FI) and the growth (or bodyweight gain) of an animal and is quantified by several indexes, such as RFI and feed conversion rate (FCR) [10][11][12].Feed efficiency is influenced by several factors, including the breed and its sex, age, diet, and management [13,14].The RFI is an important index measuring feed efficiency, and it is defined as the difference between the actual and expected feed intake [15].In 1963, the concept of RFI was first proposed in beef cattle research [16], and it was first applied to chickens by Luiting in 1991 [17], covering the calculation methods of RFI, heritability calculation, and the phenotypic and genetic correlations with relevant traits.Since RFI reflects the variation in feed efficiency, it appears to be independent of growth traits [15].Studies have shown that RFI is moderately heritable in poultry [12,18].The study by Bai et al. demonstrates that selecting for low RFI can improve poultry feed efficiency without compromising growth performance [19].Growth traits are key in poultry breeding, and the properties of growth traits need to be considered in breeding for feed efficiency.In recent years, GWAS has been applied in poultry to explore the association between host genetics and economic traits [20].Earlier studies identified candidate genes including NSUN3, EPHA6, and AGK for broilers RFI, while LAP3 was identified as a candidate gene for broiler body weight [21][22][23][24][25][26].By deeply studying the functions and regulatory mechanisms of these genes, our aim is to provide more efficient and sustainable solutions for the chicken breeding technologies [4].
In the breeding process of this study, the focus will be on maintaining the excellent meat quality while moderately improving the feed efficiency and growth rate.The objectives of this study were (1) to determine the inheritance pattern of feed efficiency and growth traits to provide a basis for formulating a better breeding method for Wenchang chickens, and (2) to find key variants affecting feed efficiency and growth traits of Wenchang chickens and also to elucidate the molecular mechanism of feed efficiency and growth traits.

Experimental birds
The study uses 1,547 chickens hatched and raised by Hainan Tanniu Wenchang Chicken Co., Ltd.These birds are from a commercial line of Wenchang chicken selected for 18 generations.All birds were raised in three-tier battery cages (one bird per cage) under the same management and nutritional conditions.The diet was formulated based on the National Research Council (NRC) requirements [27] and the Feeding Standards of Chickens established by the Ministry of Agriculture, Beijing, China [28].The chickens were raised from 87 to 113 days to collect data on individual phenotypes.

Phenotypes
The phenotypes measured in this study included body weight and feed intake of chickens.
For body weight measurements, an electronic scale with a precision of 0.1 g was utilized.The trough was emptied of any remaining feed 12 h prior to weighing, and body weights were recorded at 87, 95, 104, and 113 days of age.
Feed intake in this experiment was recorded for 26 days.Fresh feed was added at a fixed time every day and recorded.In addition to the above measured phenotypes, the phenotypes calculated from the measured phenotypes included FI, ADFI, ADG, Metabolic weight in the middle of the test (MWT), and RFI.FI, g.
RFI, g/d: The RFI was obtained from the multiple regression equation of average daily feed intake and metabolic weight in the middle of the test and average daily weight gain [18].The equation was as follows: where ADFI represents the mean average daily feed intake, µ is the intercept, sex is a fixed effect, MWT and ADG are as defined above, β1 and β2 represent the partial regression coefficient.

Genotyping, imputation, and quality control (QC)
At 112 days of age, 2 mL of blood was collected from the wing vein, placed in an anticoagulant tube containing EDTA-K2, mixed, and stored at -20℃ until later analysis.Genomic DNA was extracted from blood samples with the phenol-chloroform method.Genotyping was conducted with a customized chicken 55 Beagle 5.2 software was used to impute the 55 K chip data to the whole-genome sequence (WGS) level [30].Before imputation, inconsistencies between the target panel and the reference panel were checked using conform-gt software (http://faculty.washington.edu/browning/conform-gt.html).Then, the 55 K SNP chip data were populated to the resequencing level using Beagle 5.2 software.The QC condition: HWE at P < 1.00e-6, setting the detection rate of individual genotype at ≥ 90%, the detection rate of single SNP site at ≥ 90%, and a MAF of ≥ 0.05.A total of 12,184,765 autosomal SNPs from 1,479 samples passed QC.

Estimation of heritability and genetic correlations
Phenotypic correlation coefficients were calculated using ggpairs within the GGally R package, and then genetic parameter estimation was performed using restricted maximum likelihood (REML) of GCTA v1.93.2 beta software [31,32].The statistical model used was: where, y is a vector of observations, b is a vector of fixed effects (i.e., sex), α is the random vector representing the genomic effects, e is the vector of random residual effects, Xand Z are incidence matrices.The distribution of the random animal effect α is α ∼ N(0, Gσ 2 a ) with G is being the genomic relationship matrix, and σ 2 a being the additive genetic variance.

Genome-wide association study
A GWAS analysis was carried out between all the genotyped SNPs and feed efficiency and growth traits using a mixed linear model (MLM).The MLM for feed efficiency and growth traits was performed using sex (female or male) as a fixed effect and the top three principal components (PCs) as covariates.All association tests were performed using the MLM option in GCTA based on the following model [33]: where y is a vector of observations, X and Z are incidence matrices for the vectors for parameters b and µ, b is a vector of fixed effects including the sex and three eigenvectors from principal component analysis (PCA), µ is the vector of the additive genetic effect of the candidate SNP to be tested for association, and e is the vector of the residual effect.The Bonferroni correction method was used in this study to determine the significance thresholds, and the formula for performing the Bonferroni corrected multiple tests was as follows: Where P is the corrected significance threshold, α rep- resents the significance threshold for a single test, and N represents the number of multiple hypothesis tests, i.e., the number of SNPs analysed by GWAS.We calculated the number of genome-wide independent markers using the PLINK (V1.9) command -indep-pairwise, with a window size of 25 SNPs, a step of five SNPs, and an r 2 threshold of 0.2.Manhattan and quantile-quantile (Q-Q) plots were derived from the GWAS results using the qqman (https://cran.r-project.org/web/packages/qqman/) and Cairo (http://www.rforge.net/Cairo/)packages within R software (http://www.r-project.org/).LD blocks of target regions were performed using Haploview v4.2 software [34].For additive and dominance effects of important SNPs on traits, the calculation process in this study was done in ASReml v4.1 software (https://asreml.kb.vsni.co.uk/knowledge-base/asreml_documentation).The SNP positions were updated according to the newest release from Ensembl (https://asia.ensembl.org/index.html).Identification of the closest genes to genome-wide significant and suggestive variants was obtained using the ChIPpeakAnno package (https://www.bioconductor.org/packages/devel/bioc/vignettes/ChIPpeakAnno/inst/doc/ pipeline.html).Gene function enrichment analysis was performed using bioinformatics (https://www.bioinformatics.com.cn).

Differential expression analysis
In this study, the RFI or BW of 127 broiler chickens was ranked from low to high.Using individuals with the highest (n = 15) and lowest RFI (n = 15) phenotypes, and individuals with the highest (n = 15) and lowest (n = 15) BW.
The Wenchang chickens used in this study were a mix of males and females, so we did not differentiate between genders in the 30 chickens selected.Gene expression of important candidate genes in the different groups can demonstrate the reliability of the GWAS results.Figures were generated using GraphPad Prism 8 [4].

Description of phenotypic traits
Descriptive statistics of the feed efficiency traits of Wenchang chickens are shown in

Genetic parameters' estimates
Comparison of density markers before and after the imputation in Fig. 3 revealed a significant increase in loci after imputation.The estimated genetic parameters of residual feed intake and body weight are shown in Number of animals (N), 2 mean (M), 3 standard deviation (SD), 4 maximum (Max), 5 minimum (Min), and 6 coefficient of variation (CV) of Wenchang chickens 0.92 and 0.86, respectively (P < 0.001).Regarding phenotypic correlation, the genetic correlation between RFI and ADFI was 0.63 and decreased with body weight to a minimum of 0.20.Based on the result, the heritability of RFI in Wenchang chicken is low, while the heritability of body weight is not significantly different from that of other breeds.

Genome-wide association study of feed efficiency and growth traits
The GWAS results showed significant SNPs on GGA 2, 6 and 26 were associated with RFI (expansion coefficient λ of 0.987) (Fig. 4a-b).One of significant SNPs 6_21123592 on GGA6 was located on the intron of the PLCE1 gene, with this SNP explaining 2.46% of the genetic variation.Additionally, significant SNPs 2_45795056 and 26_2851843 were located in the lncRNA introns of the  The genetic correlation was below the diagonal and heritability of feed efficiency and growth traits were on the diagonal.The superscript ** means P < 0.01, *** means P < 0.001 ENSGALG00000052614 and ENSGALG00000001264 genes, respectively (Table 3).The GWAS results showed significant SNPs on GGA 2, 4, 6, 19, and 28 were associated with ADFI (expansion coefficient λ of 0.994) (Fig. 4c-d).Specifically, the significant SNPs 4_75971941, 4_85488222, 19_2763444, and 28_388715 were found on the genes LAP3, IMMT, GATSL2, and FBN3, respectively (Table 3).The SNP 6_21123592 for RFI in the CT genotype exhibited significantly higher values than those in the CC genotype (Table 4).Another SNP 2_45795056 showed a significant additive effect on ADFI but a nonsignificant effect on the two feed efficiency traits assessed in this study.The estimated additive effect of candidate SNP 6_21123592 on RFI and ADFI were − 0.81 and 2.67 respectively, while the dominance effect was − 2.99 and 0.93 respectively.Significant dominance effects were observed on both RFI and ADFI, while the additive effect was not significant (Table 5).
The four LD blocks were detected within the common location of growth traits on GGA4:75971.28-75974.42kb (Fig. 6), each containing 3-10 SNPs.The most significant SNP 4_75971941 was located on block1, which is located in the intron of the LAP3 gene.Genotype analysis showed that the body weight of the CC genotype at different ages was significantly higher than that of the CT genotype.Specifically, the weight of the CC genotype  The different superscripts mean the significant difference, namely P < 0.05  7).
The candidate sites for body weight were determined to be located in distinct LD blocks.In this study, haplotype association analysis was conducted using the glm model with sex as a fixed effect (Table 8).The analysis indicated that block1 had a highly significant effect (P < 0.01) on BW at various ages, while block2 exhibited a significant effect (P < 0.05) on BW across different age groups.The results revealed that block2, block3, and block4 all had a significant impact on ADFI (P < 0.05), while none of the four LD blocks showed a significant impact on RFI (P > 0.05).

Differential expression analysis
We conducted gene differential expression analysis on all genes mapped to the feed efficiency and growth traits, respectively.Significant differences in phenotypes were observed between the high and low RFI groups (Fig. 7a) and between the high and low BW groups (Fig. 7c).Between the high and low RFI groups, the PLCE1 and IMMT genes showed significantly (P < 0.01) higher expression levels in the high RFI group than in the low RFI group (Fig. 7b).Similarly, expression level of LAP3, MED28, and QDPR genes was significantly higher in the high BW group than that in the low BW group (P < 0.01), while expression of LDB2 and SEL1L3 genes was     significantly (P < 0.05) higher in the low BW group than the high BW group (Fig. 7d).

Phenotypic and estimation of the genetic parameters of Wenchang chickens
Body weight traits (BW87-BW113) of Wenchang chickens were selected for the GWAS analysis because Wenchang chickens are listing age within that age range.Through different trait measurements at this stage, our study revealed that the BW87 was 1,409.48g, that of BW113 was 1,894.36g, the ADG was 18.66 g/d, and the relative growth rate was 25.60%.These findings are consistent with previous studies on yellow-feathered chickens in China [35,36].The low heritability of RFI may be caused by the lack of systematic breeding efforts and the chickens are in the latter stage of growth, similar to the results of previous studies [37,38].Heritability of BW, ADFI and ADG revealed in the current study ranged from medium to low (0.21-0.44), which was similar to that of other breeds [39][40][41].The phenotypic and genetic correlations between RFI and ADFI were 0.92 and 0.63, similar to the research of Shirali et al [42].To optimize breeding outcomes, it is recommended to integrate additional metrics, such as ADFI, and implement a multi-trait selection approach throughout the breeding process, ultimately enhancing breeding results in Wenchang chickens.

Genome-wide association study of feed efficiency and growth traits
To investigate the genetic architecture of feed efficiency and growth traits in Wenchang chickens, a mixed linear model was used for GWAS analysis of related traits.
Through GWAS analysis, we identified significant SNP 6_21123592 on GGA6 was located within the intron of the PLCE1 gene associated with RFI.Previous studies have shown that PLCE1 gene is highly expressed in the nervous system and belongs to the phosphoinositide-specific phospholipase C family.The production of second messenger molecules such as diacylglycerol is regulated by activated phosphatidylinositol-specific phospholipase C enzymes, which mediate small GTPases of the Ras superfamily through the activity of its Ras guanine exchange factor.As the effector of heterotrimer and small G protein, PLCE1 is involved in regulating cell growth, T-cell activation, actin organization and cell survival.
Mapping the PLCE1 gene function were mostly related to nervous system activity, which regulates the function of the brain in different ways, and the brain was key to regulating diet behavior and body energy homeostasis [43][44][45][46].The significant QTL affecting BW, ADG, and ADFI was located in the 73.15-76.55Mb interval of GGA4.There were 13 important candidate genes in the colocalization interval related to BW, including LAP3, KCNIP4, NCAPG, LDB2, FAM184B, SEL1L3, ZCCHC4, PPARGC1A, PACRGL, SLIT2, LCORL, MED28 and QDPR.Among these candidates LAP3, MED28, QDPR, LDB2 and SEL1L3 demonstrated differential expression between high and low groups.Comparison with the Animal QTL database (Chicken QTL Database at (animalgenome.org)) reveals a total of 305 QTLs related to BW and 221 QTLs associated with ADG within 73.15-76.55Mb interval of GGA4.LAP3 has been shown to catalyze the hydrolysis of the amino-terminal leucine residues of protein or peptide substrates, with diverse functions in mammals, invertebrates, microbes, and plants [47].The primary function of LAP3 lies in protein maturation and degradation, processes crucial for metabolism, development, adaptation, and repair [48].LAP3 gene variation may underlie variations in growth rates among species and significant genetic polymorphism of traits of interest in breeding, potentially leading to applications in animal breeding.Furthermore, studies have been conducted on its SNP and its association with growth traits in mammals, such as bovine [49].Another study found that the LAP3 gene may have a potential function affecting muscle development in sheep [50].Prenatal development stages are directly related to the growth and development of individual skeletal muscle, which determines the number of muscle fibers and postnatal muscle mass and further exerts long-term effects on the postnatal growth of animals [51][52][53].Related research that tracked LAP3 mutations in Hu sheep populations reportedly linked, the mutations with body weight at different growth stages [54].In poultry LAP3 gene was found to be associated with chicken growth traits [22].
Related study linked the MED28 gene with live weight in sheep [55].We found that the MED28 gene was related to muscle development in pig [56].QDPR for an enzyme that regulates tetrahydrobiopterin (BH4), a cofactor for enzymes involved in neurotransmitter synthesis and blood pressure regulation.Therefore, QDPR gene are also important genes in the regulation of growth [57].Previous studies have shown that the LDB2 gene located at GGA4 is important for chicken growth traits, and a 31-bp indel was significantly correlated with multiple growth and carcass traits in the F2 population and affected the expression of the LDB2 gene in muscle tissue [26,58,59].It was also identified as an important candidate gene for rapid growth in chickens and had the strongest association with late body weight in Jinghai yellow chicken hens [22,60,61].According to relevant studies, the KCNIP4 gene is located on GGA4, and the candidate gene belongs to the potassium channel interaction protein family and has a wide range of physiological functions, including heart rate regulation, insulin secretion, neurotransmitter release, and smooth muscle contraction.It was considered to be an important candidate gene for growth traits of chickens.In addition, it was reported in different breeds and different growth stages, which also verified that the KCNIP4 and FAM184B genes can affect the growth and development of chickens [58,62,63].Studies have shown that the NCAPG gene is an important candidate gene for mammalian body size growth traits in growth trait association analysis of horses, sheep, and domestic donkeys and is involved in chromosome condensation and methylation [64][65][66][67].The FAM184B gene has been found in previous studies to be associated with cattle carcass weight [68].

Differential expression analysis
Transcriptomic data enable the quantification of DNA or RNA abundance and expression levels [69].Differential analysis between different groups is conducted by measuring the expression levels of gene RNA, thereby identifying distinct patterns and variations in gene expression among different groups.The results of related studies also confirm gene expression data of important candidate genes in different groups can demonstrate the reliability of the GWAS results [70,71].We showed that expression of PLCE1 in the high and low RFI groups of broilers, validatied it a key candidate gene for RFI.We found that the candidate gene LAP3 was related to the BW, ADFI and ADG traits, as demonstrated by the significant difference in the expression of LAP3 in high and low body recombination in broilers, and relevant research reports also indicated that this might be an important candidate gene affecting growth traits.This study also found that the candidate genes LDB2, KCNIP4, FAM184B, and NCAPG were located on GGA4 and were related to growth traits, which provides an important reference value for subsequent research on the growth traits of Wenchang chickens.
The validity of the SNPs and candidate genes obtained in this study is worth further extensive verification.The significant SNPs and candidate genes identified in this study can be incorporated into the chip markers in future research.By utilizing genomic selection breeding techniques, these markers can be used to breed and improve the target traits of WenChang chickens, ultimately enhancing breeding efficiency.

Conclusions
In this study, we identified that the significant SNP 6_21123592 was located in candidate gene PLCE1 for feed efficiency traits of Wenchang chickens, and the significant SNP 4_75971941 was located in candidate gene LAP3.Other candidate genes including MED28, QDPR, LDB2, and SEL1L3 were identified for growth traits of Wenchang chickens.This provides a good theoretical basis for developing methods of Wenchang chickens breeding, and by further studying the functions and regulatory mechanisms of these genes, we could provide more efficient solutions for breeding of these chickens.

F
I = f eed amount duringthe test − W eight of f eed remaining at the end of the test (1) ADFI, g/d.

Fig. 3 Fig. 2
Fig. 3 Comparison of density markers before and after the imputation

Fig. 4
Fig. 4 Manhattan and Q-Q plots of GWAS for feed efficiency traits.The horizontal red lines indicate the thresholds for genome-wide significance (P = 3.47E-08), and the horizontal blue lines indicate the thresholds for suggestive significance (P = 6.93E-07)

Fig. 5
Fig. 5 Manhattan and Q-Q plots of GWAS for growth traits.The horizontal red lines indicate the thresholds for genome-wide significance (P = 3.47E-08), and the horizontal blue lines indicate the thresholds for suggestive significance (P = 6.93E-07)

Fig. 7
Fig. 7 The expressions of candidate genes for feed efficiency and growth traits.(a), Analysis of RFI Differences between high and low RFI groups.(b), Differential expression of related candidate genes between high and low RFI groups.(c), Analysis of BW Differences between high and low BW groups.(d), Differential expression of related candidate genes between high and low weight groups.Data were expressed as the mean ± SEM (n = 15), *** P < 0.001, ** P < 0.01, * P < 0.05

Table 2
. The heritability of the RFI, BW, ADFI and ADG ranged from moderate and low level (0.05-0.44).RFI had the highest genetic association with ADFI and ADG, at

Table 1
Descriptive statistics for feed efficiency and growth traits of Wenchang chickens (N 1 = 1532)

Table 2
Estimates of genetic parameters for feed efficiency and growth traits of Wenchang chickens

Table 3
The significant SNPs information of the feed efficiency traits in Wenchang chicken

Table 4
Effect of candidate SNPs on feed efficiency traits

Table 5
Additive and dominance effects of candidate SNPs on feed efficiency traits1

Table 6
The significant SNP information of growth traits in Wenchang chicken

Table 7
Additive and dominance effects of candidate SNPs on ADFI and growth traits The superscript * means P < 0.05, ** means P < 0.01

Table 8
The effect of different LD blocks on feed efficiency and growth traits in Wenchang chickens A significant impact on traits is indicated by P < 0.05 for LD blocks